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in ; ABSTRACT 

! We report the results of ID particle-in-cell simulations of ultrarelativistic 

shock waves in proton-electron-positron plasmas. We consider magnetized shock 

. waves, in which the upstream medium carries a large scale magnetic field, directed 

transverse to the flow. Relativistic cyclotron instability of each species as the 
incoming particles encounter the increasing magnetic field within the shock front 
provides the basic plasma heating mechanism. The most significant new results 
come from simulations with mass ratio m p /m± = 100. We show that if the 
protons provide a sufficiently large fraction of the upstream flow energy density 
(including particle kinetic energy and Poynting flux), a substantial fraction of 
the shock heating goes into the formation of suprathermal power-law spectra of 
e~ — e + . Cyclotron absorption by the pairs of the high harmonic ion cyclotron 
. J_h ! waves, emitted by the protons, provides the non-thermal acceleration mechanism. 

As the proton fraction increases, the non-thermal efficiency increases and the 
e~ — e + power-law spectra harden. 

When the proton fraction is small (pair plasma almost charge symmetric), the 
e~ and e + have approximately equal amounts of non-thermal heating. At the 
lower range of our simulations with mass ratio 100, when the ions contribute 56% 
of the upstream flow energy flux, the pairs' non-thermal acceleration efficiency 
by energy is about 1%, increasing to 5% as the ions' fraction of upstream flow 
energy increases to 72%. When the fraction of upstream flow energy in the ions 
rises to 84%, the efficiency of non-thermal acceleration of the pairs reaches 30%, 
with the e + receiving most of the non-thermal power. 

We suggest that the varying power law spectra observed in synchrotron 
sources that may be powered by magnetized winds and jets might reflect the 
correlation of the proton to pair content enforced by the underlying electrody- 
namics of these sources' outflows, and that the observed correlation between the 
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X-ray spectra of rotation powered pulsars with the X-ray spectra of their nebulae 
might reflect the same correlation. 

Subject headings: acceleration of particles — shock waves — ISM:jets and out- 
flows 



1. Introduction 

Collisionless shock waves in relativistic astrophysical flows have been implicated as the 
energizing mechanism for emission in a variety of non-thermal sources - termination shocks of 
pulsar winds (Slane (2005); Kaspi et al. (2004)), hot-spot shocks terminating jets from radio 
galaxies (Heavens & Meisenheimer (1987); Wilson, Young & Shopbell (2000); Balucirinska- 
Church et al. (2005)) and internal shocks in AGN and microquasar jets (Marscher & Gear 
(1985); Tiirler et al. (2000, 2004)), internal and external shocks in GRBs (Piran (2004)) 
and outflows from other high energy transients (e.g., Gaensler et al. (2005)). While shock 
jump conditions applied to models of magnetohydrodynamic flow set the constraints on 
what fraction of flow energy goes into heat of all sorts, the partition of that energy between 
ions and electrons (and positrons, when these are present), and the distributions of particle 
momenta (thermal or non-thermal), are kinetic questions not addressed by the macroscopic 
conservation laws. The kinetic questions are those of collisionless plasma physics - two body 
encounters cannot generate the entropy implied by the jump conditions on the observed 
scales of transition from non-radiative flow (upstream, in the shock excitation model) to 
the observed regions of synchrotron emission (downstream, in the shock model), nor can 
collisional processes lead to non-thermal particle spectra. 

A fundamental approach in which the flow and kinetic structure of relativistic shocks 
are resolved on all scales offers the opportunity to evaluate the nature of the electromagnetic 
turbulence generated, as a function of the upstream parameters of the flow, simultaneously 
with a determination of the heated (possibly non-thermally accelerated) particle spectra 
both up and downstream. For this purpose, fully kinetic particle-in-cell (PIC) simulations 
(Birdsall & Langdon (1991); Dawson (1995)) provide a powerful investigatory tool, since 
these model the plasma from first principles, with the only approximations being those that 
go into the method - primarily the "cloud-in-cell" algorithm, which in effect softens the 
interactions to the grid scale. 

PIC simulations of relativistic transverse magnetosonic shocks in a symmetric pair 
plasma were first reported by Langdon et al. (1988) and by Gallant et al. (1992), using 
spatially ID models. Magnetized shocks are characterized by coherent gyration (Alsop & 
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Arons (1988)) and cyclotron instability (Hoshino & Arons (1991)) of the plasma species 
as they encounter the shock front - magnetic reflection mediates the density and velocity 
transition, while the instabilty of the induced gyration mediates the entropy generation. Sim- 
ulations of the Weibel instability when B = with applications to unmagnetized relativistic 
shocks have been reported recently by a number of authors (Silva et al. (2003); Frederik- 
sen et al. (2004); Hededahl et al. (2004); Jaroschek (2005); Hededahl k Nishikawa (2005)). 
Spitkovsky and Arons (in preparation) have found comparable results, and show how the 
shock structures in pair plasmas depend on the magnetic field strength and orientation with 
respect to the flow. All these results indicate that non-thermal particle acceleration is either 
weak or absent in shocks in pair plasmas. 

Hoshino et al. (1992) found that when the upstream flow includes heavy ions which 
carry a significant fraction of the upstream flow energy, non-thermal acceleration does occur, 
as a consequence of the absorption of ion waves emitted at harmonics of the ion cyclotron 
frequency reabsorbed by the pairs, even though the shock structure is one dimensional and 
no cross field diffusion occurred. The interaction is resonant, in contrast to the non-resonant 
scattering underlying Diffusive Shock Acceleration (DSA). Computational limitations re- 
quired the use of a low rest mass ratio, mi/m± = 10, 20. 

Those inhomogeneous shock calculations showed no electron acceleration. While most 
of the characteristics of the instability scale with just the energy density ratio between the 
protons and the pairs, the efficacy of the waves as accelerators of electrons and positrons 
also depends on the polarization properties of the proton cyclotron waves. The lack of 
electron acceleration at low mass ratio stems from the polarization of the waves excited. In 
a quasi-neutral pair plasma the waves are linearly polarized. If the heavy ions are a small 
fraction by number, the waves excited by the ion maser remain almost linearly polarized. 
Since linearly polarized waves are an equal mixture of left and right circularly polarized 
modes, equal electron and positron acceleration might have been expected. However, the 
simulations also revealed that the mechanism is a significant accelerator only if the energy 
density in the ions Uu = uiiNn^i — l)c 2 , where Nn is the upstream ion density and 71 is 
the upstream flow Lorentz factor, is comparable to the upstream energy density in the pairs, 
U 1± = m±( 7 i - l)c 2 (N 1+ + Aq_). Then N U /(N 1+ + N,_) = (mJm^Uu/U^) ~ m ± /m,. 
For the small mass ratios employed, the dispersion relations for small amplitude waves show 
that at comparable energy densities in the species, the waves are still strongly left handed, 
with consequent preference for positron acceleration. 

The earlier work also entirely neglected upstream thermal spread in the momenta of the 
species. The acceleration mechanism relies upon cyclotron absorption of waves in the pairs 
at high harmonics of the ion cyclotron frequency. Since the shock in the pairs alone heats 
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the pairs to a temperature T 2 ± ps m±c 2 7i (if <7i = -B 2 /47rpi7ic 2 <C 1, with £?i the upstream 
magnetic field strength and p 1 the upstream rest mass density, both measured in the down- 
stream frame), efficient non-thermal acceleration of the pairs requires cyclotron absorption 
of ion waves at frequency uo m Q c ± = eB /m±c~fi = (m i /m±)eB/m i c'~fi = (m,i/m±)VL ci \ that 
is, the relativistic cyclotron instability of the ion ring in the shock front's leading edge must 
generate significant power at harmonic numbers n ~ rrii/m± ^> 1. When the upstream 
plasma is completely cold, this is not a problem - the growth rate of the instability in a 
homogeneous medium scales as Vl ci n~ 1 / 9 (Hoshino & Arons (1991)). We shall see that if 
there is thermal dispersion in the upstream ion momenta Ap, growth is inhibited, with no 
substantial wave power at harmonic numbers above rriiCji/ Ap. 

In this paper we reinvestigate the one dimensional shock structure of a transverse mag- 
netosonic shock, with particular attention to the effect of more realistic mass ratios and the 
incorporation of thermal dispersion in the upstream medium. 

In §2 we undertake the study of the structure and acceleration properties of transverse 
relativistic shocks. §2.1 contains a study of shocks in pure pair plasma. Consistent with 
previous results, we find no signs of non-thermal particle acceleration. We introduce a 
proton component in §2.2 and find, for the first time in the literature, evidence of non- 
thermal acceleration not only of positrons but also of electrons. In §2.3 we show that at high 
mass ratio the waves become close to linearly polarized when the ion and pair energy densities 
are comparable. This leads to comparable heating and acceleration rates for positrons and 
electrons. In §2.4 we include in the shock simulations a thermal spread of the ions' upstream 
distribution and study how this affects the results. Finally, we present a summary of our 
results and our conclusions in §3. Aspects of the linear stability of relativistic protons 
gyrating in a relativstically hot pair plasma are described in the Appendix, §A.l. 

2. Structure and Non-thermal Particle Acceleration Properties of Relativistic, 

Transverse Magnetosonic Shocks 

We study relativistic transverse magnetosonic shock waves. Our analysis was carried out 
through numerical simulations performed with XOOPIC (Verbonocoeur, Langdon & Gladd 
(1995)). The simulations were ID and were performed in the geometry shown in Fig. 1 
below. 

When the simulation starts, the simulation box is permeated with a static uniform 
magnetic field directed along the z-axis and a uniform electric field along the y-axis. New 
plasma is constantly injected in the simulation box from the left boundary. The plasma 
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moves along the x-axis with velocity v xl (four- velocity u x i). The drift motion is initialized 
and maintained in the upstream medium by E y , which is set equal to E = — v xl /c x B. 

The right boundary of the simulation box is a conductor which behaves as a perfectly 
reflecting wall, both for particles and electromagnetic radiation. When the plasma first 
impacts on the wall, a reverse shock propagates towards the left. What we are interested in 
is the behaviour of the fluid at, and right after, the crossing of this reverse shock. 



shock front 




upstream 



downstream 



Fig. 1. — The simulation geometry. 



The upstream and downstream conditions are related by the ideal MHD jump condi- 
tions for a relativistic, transverse, magnetosonic shock. These were derived by Kennel & 
Coroniti (1984), as expressed in the shock frame, and by Gallant et al. (1992) as expressed 
in the downstream frame. Quantities expressed in this latter frame are more convenient for 
purposes of direct comparison of the results of our numerical simulations with the theoretical 
expectations. Since the wall on the right of the simulation box is at rest, the ^-component 
of the fluid velocity must vanish downstream of the shock. The downstream frame is there- 
fore coincident with the computational observer's frame and it is the frame in which all 
quantities that will enter the following discussion are measured. In this frame, the shock 
moves toward the left, with velocity —cj3 s hock- We assume the fluid upstream to be cold, 
Pi riiinc 2 , where pi and n\ are the proper pressure and density upstream, and to move 
at ultrarelativistic speed, with Lorentz factor ji ^> 1 and 4- velocity cu\ ~ c^±. For the fluid 
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downstream, we assume the ideal ultrarelativistic equation of state to hold, with proper 
enthalpy w 2 ~ Tp 2 /(T — 1). The adiabatic index T would be equal to 4/3 for an isotropic 
relativistic fluid. But in our analysis the particles' momenta are restricted to the plane 
perpendicular to B therefore the appropriate value of T is 3/2, describing a 2-dimensional 
ultrarelativistic fluid. 

After expressing the densities in terms of the simulation frame values N = wy, the jump 
conditions reduce to: 
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where we have used the definition of the magnetization parameter 
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Using Eq. 3 and Eq. 4 to eliminate p 2 and B 2 from Eq. 2, one finds 
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The jump in all quantities at the crossing of the shock front is then determined as a function 
of a and T. The density and magnetic field compression are determined through Eqs. 1 and 
4, while the ratio between the thermal pressure downstream and the ram pressure upstream 
can be similarly expressed in terms of a and /9 s hock after manipulation of Eq. 2. One obtains: 
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(6) 



The low and high a limits (ct <^ 1 and a ^> 1, respectively) of Eqs. 1, 4, 5 and 6 are 
discussed extensively by Gallant et al. (1992), and the main result is that with increasing 
magnetization the fraction of upstream energy that goes into an increase of the magnetic 
pressure downstream becomes larger, at the expense of the particle thermal energy. In the 
following, we shall deal with the case of a plasma with a ~ 1. 
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2.1. Numerical simulations of relativistic transverse shocks in a pure pair 

plasma 

The results we show and discuss in the following refer to a simulation performed on a 
plasma with a_ = a + = 2, corresponding to a tot = 1, and Lorentz factor upstream = 40. 
The simulation grid was made of 1024 cells, each of size Ax = r L /10, with r L being the 
particle Larmor radius based on the upstream B field and flow Lorentz factor. The initial 
number of particles per cell was 32 for each species. The results are basically unchanged 
when Ax is varied from r^/10 to the number of cells between 512 and 2048, and the 

number of particles of each species per cell from 8 to 64. The time-step is At = 0.9Ax/c in 
all cases, in order to satisfy the Courant condition. 




x/r L 



x/r L 



x/r L 



x/r L 



Fig. 2. — Plots on the left: The electron four-velocity as a function of position in the 
box. The x-coordinate is in units of t\l, the pair Larmor radius upstream. The shock front 
is located at x s /r^ w 50. Top panel: x-component of the four velocity. Middle panel: y- 
component of the four-velocity. Bottom panel: a zoom of the top panel onto the region 
around the shock front. Plots on the right: The electromagnetic field as a function of 
position in the box. The x-coordinate is again in units of the pair Larmor radius upstream. 
Top panel: longitudinal component of the electric field. Middle panel: transverse component 
of the electric field. Bottom panel: the magnetic field. 

As it was the case for the simulations performed by Gallant et al. (1992), we find in 
our simulation the presence of an electromagnetic precursor wave, propagating ahead of 
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the shock, into the upstream medium, at essentially the speed of light. The results in the 
following figures refer to a time when the reverse shock has propagated far enough from the 
right wall of the simulation, so that edge effects (spurious oscillations due to the vicinity 
of the conducting boundary) do not affect the dynamics around the shock front, while the 
precursor has not reached the left boundary of the box yet. 



40<x/r L <45 



40<x/r L <45 





Fig. 3. — The electron (panels on the left) and positron (panels on the right) distribution in 
u x — u y space in different regions around the shock front, as specified on top of each panel. 
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In Fig. 2 we show, on the left, the (u x — x) and (u y — x) sections of the electrons' 
phasespace. The behavior of the positrons is completely analogous. The main thing to 
notice is the loop that develops at the crossing of the shock in the (u x — x) plane, especially 
evident in the bottom panel. This is the initial reflection loop described by the soliton 
solution of Alsop & Arons (1988). 

The particles' four-velocities in the reflection region form distorted cold rings in mo- 
mentum space, as shown in Fig. 3, so that one can expect cyclotron instability to develop, 
as expected from the uniform medium instability theory outlined in the Appendix with the 
ions omitted (see also Hoshino & Arons (1991)). 

In Fig. 2 we also plot the different electromagnetic field components as a function of 
position. The longitudinal electric field always stays around zero with negligibly small oscil- 
lations. This is just what one expects due to the charge and mass symmetry, together with 
overall neutrality, leading to decoupling of the transverse from longitudinal extraordinary 
mode. 

The transverse component of the electric field vanishes on average behind the shock, 
as predicted by ideal MHD, but of course there are oscillations due to the development of 
cyclotron instability in the pairs. The same behavior is seen in the magnetic field. The latter, 
after the first two main peaks, corresponding to the turning points of the first reflection loop, 
oscillates around the value computed from the jump conditions (Eqs. 4 and 5) for a plasma 
with (7 = 1 and T = 3/2 (dashed curve in the plot). 

Finally, in Fig. 4 we show the evolution of the pair distribution function at the crossing 
of the shock. It is apparent that no signs of non-thermal acceleration are observed. The 
distribution function progressively evolves towards a maxwellian for both species. The tem- 
perature of the maxwellian is in perfect agreement with the prediction of ideal MHD (Eq. 6): 
T eff = kT 2 /m±c 2 « 0.4 7l . 
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Fig. 4. — The evolution of the particle distribution after the crossing of the shock. The 
different line-types refer to particles extracted from different slices of the simulation domain 
as specified in the figure: short-dashed for particles at positions 40 < x/tl < 50, long-dashed 
for positions 50 < x/r L < 60 and solid for positions 60 < x/r L < 80. The top panel refers to 
electrons and the bottom one to positrons. In both cases 7(7) is normalized to the particle 
number denisty upstream. The thick solid curve represents a 2-D relativistic maxwellian of 
temperature T e g (see text). 
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2.2. Numerical Simulations of relativistic transverse shocks in e — e + — p 

plasmas: acceleration of electrons 

We show that with protons present in the upstream flow, the downstream exhibit 
non-thermal heating, a.k.a. shock acceleration. We first tackle the problem of how the 
particle acceleration process is affected by the waves' polarization. 

In the following we discuss the results of two simulations in which a plasma with a_ = 2 
and <7j = 0.25 is considered. The latter condition is realized in two different ways, first 
adopting a mass and number ratio between ions and electrons equal to 20 and 0.4 respectively 
(case A), and then using nii/m± = 40 and iVj/AL = 0.2 (case B). 

The simulation geometry is that in Fig. 1. In case A the simulation grid contains 2048 
cells, each of size rn±/10, where tl\± = m±c 2 ^i/ eB\, while in case B the cell size is the 
same but the number of cells is doubled. In both cases the box contains more than 10 ion 
Larmor orbits and the number of particles per cell is 16. The initial plasma Lorentz factor 
is 71 = 40, but no dependence of the results on this parameter is expected as long as the 
condition 71 3> 1 is satisfied. As to the time-step, finally, this is taken to be At = 0.9Ax/c, 
so that the Courant condition is ensured. 

In Figs. 5-6 we show, on the left, the section (u x — x) of the phase-space of the different 
species in the two different cases. The reverse shock is at around x s « 100r L in Fig. 5 and 
x s 200r L in Fig. 6. 

Comparing the left panels of the two figures, the first ion reflection loop is well defined 
in both cases, and the overall velocity profiles look very similar. At a closer look, however, 
one notices that in case B (Fig. 6) the spread in the electron four- velocity post-shock is larger 
than in case A, while the opposite happens for the positron velocities. This is a first signature 
of what we will find to be the most important difference between the two simulations. 

Before discussing this point, however, let us briefly comment on the plots showing the 
electromagnetic field components in the two cases. First of all, in both cases, the presence 
of the electromagnetic precursor propagating ahead of the shock in the upstream medium 
is well evident. This wave is a purely transverse electromagnetic wave, as one can readily 
deduce from the absence of fluctuations in the longitudinal component of the electric field 
upstream of the shock. It is the same precursor as appears in magnetized shocks in pure 
pair plasmas (see previous section and also Gallant et al. (1992)), and is emitted by the 
quasi- coherent rings of e ± reflected at the leading edge of the whole shock structure. 

As to the fields downstream, now the oscillations observed in the longitudinal component 
of the electric field are larger than in Fig. 2: due to the fact that the charge and mass 
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Fig. 5. — Left panel: The (u x — x) section of the phase space of the different species for case 
A: electrons in the top panel, positrons in the middle panel and ions in the bottom panel. 
Right panel: The electromagnetic field components in the simulation labeled A: the fields 
are normalized to the upstream magnetic field intensity. In all plots the x-coordinate is in 
units of the pair gyration radius upstream. 
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Fig. 6. — Same as Fig. 5 but for case B. 
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Fig. 7. — The downstream power spectra of the different electromagnetic field components, 
for case A (dashed curves) and case B (solid curves). The power spectrum is normalized to 
the value of the magnetic field upstream, while the wavenumbers are in units of fi ce /c, with 
Vt ce being the pairs' Larmor frequency upstream. 
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symmetry of the pure pair plasma is broken by the presence of the ions, the transverse mode 
is partly longitudinal. The magnetic field oscillates around the enhanced post-shock value 
predicted by ideal MHD. The properties of the fields downstream are similar in the two cases 
we are now considering as it is clear from Fig. 7, where we compare the downstream power 
spectra extracted from the two simulations. 

In both cases one finds that the power as a function of frequency can be approximated 
around the pair cyclotron frequency as a power law with index -2. 




Fig. 8. — The postshock evolution of the distribution function of the different species, each 
normalized to the corresponding upstream density. The left panels refer to case A and the 
right panels to case B. The different line-types refer to different ranges of particle positions, 
as speciefied in the lower panel of each column. For electrons and positrons (upper and 
middle panels of each column) two fits of the distribution functions are presented together 
with the data. The thick dashed curve is a 2-D relativistic maxwellian with temperature T", 
while the thick solid curve is a maxwellian plus power- law fit (see text). 

Finally we discuss the postshock evolution of the particle distribution functions. In 
Fig. 8 we plot for both cases (case A on the left and case B on the right) the distribution 
function of the different species within different regions behind the shock. 

We observe, for all species, and in all cases, a progressive broadening of the distribution 
function with increasing distance from the shock front. However the details of the process 
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are different between the different species and between the two simulations. 



In case A, the energy release of the ions is quicker than in case B. One also sees, in the 
plots on the left, that a preferential heating of positrons occurs, as shown by the effective 
temperatures of the different species (value of T eS reported in the plots). Moreover, positrons 
are not only heated but also non-thermally accelerated, since the energy excess is found in 
the form of a well developed high energy suprathermal tail. 

In case B, one still observes hotter positrons than electrons, but the difference between 
the effective temperatures of the two species is smaller. The final value of T eff is, for both 
species, smaller than in case A, and the suprathermal positron tail is reduced, despite being 
still evident, while signs of electron acceleration are found for the first time. 

We fit the final distribution functions of the pairs (thick solid curves in the top and 
middle panels of Fig. 8) with a function N("f) given by the superposition of a 2-D relativistic 
maxwellian plus a power-law tail, with the latter exponentially cut-off at low energies: 
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In Eq. 7, N\ is the upstream density of the species considered and T' = T/m±c 2 , a and r are 
the power-law indices of the suprathermal tail and of the cut-off, respectively, and, finally, 
7 m is the Lorentz factor at which the power-law tail is cut-off. 

The best fit function is plotted in each of the relevant panels of Fig. 8 as a thick solid 
curve, while the thick dashed curve represents a maxwellian of temperature T' = T/m±c 2 . 
The curves we plot are defined by using for the parameters appearing in Eq. 7 the values 
listed in Table 1. In the last row of the same table we report the acceleration efficiency in 
terms of energy, r\. This is computed as the ratio between the energy density of e~ and e + 
in the suprathermal tail of the distribution, U nt ± and the total upstream flow energy density 

v± = u nt ±/u tot , 
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and U nt ± computed as the integral of the second term in Eq. 7 multiplied by m±7c 2 : 
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In the latter equation we see that a parameter ^ ma x appears, that refers to the maximum 
Lorentz factor up to which particles are accelerated. The theoretical upper limit for 7 max 
within the framework of the proposed acceleration mechanism is 7^ tt3 , = m i /m±'ji, i.e. the 
pairs cannot be accelerated up to energies exceeding the ions' energy upstream. 



Simulation parameters 




case A 


case B 


m,i/m± 


20 


40 


Ni/N- 


0.4 


0.2 


Ui/Utot 


0.69 


0.68 


Ui/U- 


8 


8 


Downstream e ± distribution parameters 




electrons 


positrons 


electrons 


positrons 


f 


1 


0.6 


0.92 


0.8 


T 


25 


25 


17 


16 


1m 




100 


80 


80 


N t 





3 x 10~ 3 


1 x 10~ 3 


2 x 10" 3 


a 




1.7 


2.2 


1.8 


r 




4 


4 


4 


Imaxl Imax 




1 


0.3 


0.75 


V [%} 


< 2 


20 


3 


11 



Table 1: In the upper part of the table we report the values of the paramaters with which the 
simulation was performed: mass ratio between the species; number ratio between ions and 
electrons; fraction of the total upstream energy density carried by the ions; ratio between the 
ions' and electrons' energy density upstream. In the lower part of the table, we report the 
best fit values of the parameters in Eq. 7, describing the electrons' and positrons' distribution 
downstream. In the bottom line we report the energy acceleration efficiency for the pairs, 
estimated (Eq. 8) as the percentage of the total upstream energy density U to t converted into 
non-thermal e ± in each case. 

The results presented in this section are summarized as follows. When the ions are only 
20 times more massive than the pairs and 40% of the total number of positive charges, only 
positron acceleration is observed. When a different simulation plasma is considered, with 
the same amount of energy in baryons as before, but given to a smaller fraction by number 
of more massive ions, the fraction of energy that is converted into accelerated positrons is 
found to be smaller; their spectrum is steeper and also the maximum Lorentz factor reached 
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is lower. At the same time, however, signs of electron acceleration start appearing, although 
the efficiency of this latter process is lower compared to that of positron acceleration; the 
e~ spectrum is softer and their maximum Lorentz factor smaller than that of e + . It is 
important to notice that while the results concerning simulation A are in agreement with 
those of Hoshino et al. (1992), the simulation B is the first of this kind in the literature to 
show traces of electron acceleration, thanks to the larger mass ratio and correspondingly 
lower number ratio between the ions and the pairs employed. This confirms the role of the 
wave polarization in the acceleration process. 



2.3. Relativistic transverse shocks in e — e + — p plasmas: rrii/m± = 100 

We now analyze how the situation changes increasing the mass-ratio between the ions 
and the pairs to mi/m± = 100. 

In Fig. 10 we show the energy distribution of the particles at different positions behind 
the shock, for three simulations with rrii/m± = 100 and iVj/AL varying between 0.05 and 
0.2. In all simulations we employed a box of 10240 cells and the number of particles of each 
species per cell was 12. The size of each cell was Ax = rj, e /10, with r^ e being the pair 
Larmor radius upstream. The time-step was At = 0.9 Ax /c. 

The snapshot of the simulation to which the plots in Fig. 10 refer was taken when the 
shock position was around x s /rL e = 500. The particle distribution in velocity space and 
the electromagnetic field components are shown in Fig. 9 for the case iVi/AL = 0.1. From 
comparison of the left panel of this figure, with the left panels of Figs. 5 and 6 one sees that 
the electrons show an increasingly large spread in their velocity distribution with decreasing 
Ni/N_. This feature is made more quantitative by Fig. 10 and Table 2, where simulations 
performed with the same value of the mass ratio (wij/m± = 100) but different fractions of 
ions are compared. 

It is apparent from Fig. 10 that efficient acceleration of both electrons and positrons is 
found for Ni/N- > 0.1. A very large acceleration efficiency is observed in the case when 
Ni/N_ = 0.2, for both species of particles: 50% of the positrons and 35% of the electrons are 
found in the suprathermal tail of the distribution in this case, with an energy acceleration 
efficiency that almost reaches 30% for positrons. The ratio between the efficiencies of electron 
and positron acceleration is found to be similar to that in simulation B (m,/m ± = 40, 
Ni/N- = 0.2) of the previous section, so that the conclusion we derive is that this is strictly 
related to the degree of neutrality of the pair plasma, which affects the polarization of the 
waves. 
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x/ r L x/ r L 

Fig. 9. — Left panel: The (u x — x) section of the phase space of the different species for the 
simulation with mi/m± = 100 and iVj/iV_ = 0.1. Right panel: The electromagnetic field 
components, normalized to the upstream magnetic field intensity, for the same simulation. 
In all plots the ^-coordinate is in units of the pair gyration radius upstream. 




Fig. 10. — The postshock evolution of the distribution function of the different species, 
each normalized to the relative upstream density. Left column : Ni/N_ = 0.05. Middle 
column: iVj/AL = 0.1. Right column: Ni/N- = 0.2. The notation is the same as for 
Fig. 8, with the thick solid curve representing a fit of the distribution with the function 
described by Eq. 7. The best fit values of the parameters are those in Table 2. 
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Simulation parameters 


rrii/m± = 100 


N t /N_ 


0.05 


0.1 


0.2 


Ui/Utot 


0.56 


0.72 


0.84 


Ui/U_ 


5 


10 


20 


Downstream distribution parameters 




electrons 


positrons 


electrons 


positrons 


electrons 


positrons 


f 


0.97 


0.97 


0.7 


0.7 


0.65 


0.5 


T 


18 


16 


18 


17 


15 


18 




100 


100 


100 


90 


100 


85 


N t 


3.5 x 10" 4 


3.7 x 10" 4 


2.2 x 10~ 3 


2.3 x 10" 3 


3.2 x 10~ 3 


4.2 x 10" 3 


a 


3.5 


3.5 


3.2 


2.9 


2.7 


1.6 


r 


4 


4 


3 


3 


3 


3 


Umax/ Imax 


0.15 


0.3 


0.2 


0.4 


0.3 


0.8 


V [%] 


0.9 


1 


5 


5 


5 


27 



Table 2: The table refers to simulations employing a mass-ratio of 100 between the ions and 
the pairs and number ratios N p i/N e i varying between 0.05 and 0.2. In the upper part of 
the table we report the values of the paramaters with which each simulation was performed: 
number ratio between ions and electrons; fraction of the total upstream energy density 
carried by the ions; ratio between the ions' and electrons' energy density upstream. In the 
lower part of the table, we report the best fit values of the parameters in Eq. 7, describing 
the electrons' and positrons' distribution downstream. In the bottom line we report the 
energy acceleration efficiency for the pairs, estimated (Eq. 8) as the percentage of the total 
upstream energy density U to t converted into non-thermal e ± in each case. 
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The percentage of energy in suprathermal particles is drastically reduced when the 
number of ions is reduced from iVj = 0.2A_ to iVj = 0.1 iV_ and as a consequence the energy 
ratio is halved. At the same time, the pair plasma becomes closer to neutral and the behavior 
of the different species of light particles becomes increasingly similar. For iVj = 0.1 N_ the 
energy content of the high energy tails of the distributions of the two species is the same, 
and the slopes of the suprathermal power laws become similar, though steeper than when 
the ions are twice as many. 

If Ni/N_ is decreased further by a factor of two the acceleration efficiency is reduced by 
a factor of 5 for both species and the power law indices increase further, becoming exactly 
the same. 



2.4. Relativistic transverse shocks in e — e + — p plasmas: Effect of Finite 



The results described in §2.2 and §2.3 assumed the upstream flow to be completely 
cold, i.e., to have zero momentum dispersion. The linear theory outlined in §A.l suggests 
that if the upstream flow is warm enough, magnetosonic waves at high harmonics of the 
ion cyclotron frequency will be absent in the medium behind the shock in the pairs, thus 
reducing the efficiency of injection of pairs into the acceleration mechanism. The lack of 
the harmonics above n = Ui/Au, where Au is the momentum dispersion of the upstream 
ion flow, means that the minimum energy that can be accelerated by cyclotron resonance 
is m±c 2 'y m i n = (rriiC 2 / 'M s )ji, if M s = u±/ Au < m,i/m±. Thus the suprathermal tail starts 
farther and farther out on the Maxwellian tail of the pairs, for smaller and smaller M s , thus 
reducing the acceleration efficiency (and steepening the power law part of the spectrum). 
The maximum energy possible is still set by resonance with the fundamental, but now there 
are many fewer such particles. 

The simulations whose results are shown in the following were performed assuming 
gaussian functions to describe the particles' distribution upstream of the shock: 



M s = -ui /Am is the relativistic flow sonic Mach number as defined in the linear instability 
theory. 

All three species of particles were injected in the box according to the distribution 
function in Eq. 11, with u\ corresponding to a Lorentz factor 71 = 40 and Au = 0.1 u\. The 
mass ratio between the ions and the pairs was taken to be m,i/m± = 100 and two different 



Upstream Temperature 




(11) 



- 21 - 



simulations were performed, corresponding to Ni/N^ = 0.1 and iVj/AL = 0.2. Except for the 
substitution of a Dirac 8 with a gaussian as the function describing the particles distribution, 
the simulation set-up is exactly the same as that described in §§2.3, i.e. same number of cells 
and of particles per cell, same cell-size (in units of the mean pair Larmor radius upstream) 
and time- step. 

Fig. 11 is a snap-shot of the simulation with iVj/AL = 0.1, taken when the shock front 
is at x/tl — 500. In this figure we plot the particles' phase-space and the electromagnetic 
fields in the simulation for a direct comparison with those in Fig. 11, which referred to the 
same incoming plasma taken to be completely cold. 






- — ^\/%r^- 







500 



500 

x/r. 




500 

x/r L 



x/r L 




x/r L 



Fig. 11. — Left panel: The (u x — x) section of the phase space of the different species for a 
simulation with rrii/m± = 100 and Ni/N_ = 0.1. The particle distribution function upstream 
of the shock is described by Eq. 11 with u\ ~ 40 and Au = 0.1 U\ (M s = Ui/Au = 10). 
The shock front is at x/tl — 500. Right panel: The electromagnetic field components, 
normalized to the upstream magnetic field intensity, for the same simulation. In all plots the 
re-coordinate is in units of the mean (corresponding to 71) pair gyration radius upstream. 



From the comparison of Fig. 11 with Fig. 9 it is immediately apparent that the maximum 
pair energy downstream is reduced when a spread in the particle velocities upstream is 
introduced. 

The differences are even more evident from the comparison between Fig. 10 and Fig. 12. 
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In the latter we plot the particle distribution functions as extracted from slices of the sim- 
ulation at different distances from the shock front. The left and right panels refer to the 
case when Ni/N_ = 0.1 and Ni/N_ = 0.2 respectively, and hence should be compared to 
the middle and right panel in Fig. 10. The notation is the same as for that figure. 

The particle acceleration efficiency is reduced with respect to the case when the upstream 
plasma was taken to be completely cold. As one can see by comparing the values of the 
parameters in Table 3 with those in Table 2, for both values of iVj/7V_, the fraction of energy 
carried by suprathermal particles is smaller and the power-law index of the suprathermal 
tail is steeper. 
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Fig. 12. — The postshock evolution of the distribution function of the different species, 
each normalized to the relative upstream density. Left column : iVj/AL = 0.1. Right 
column: Ni/N_ = 0.2. The notation is the same as for Fig. 10, with the thick solid curve 
representing a fit of the distribution with the function described by Eq. 7. The best fit values 
of the parameters are those in Table 3. 



The resonance condition for absorption of an ion cyclotron wave by an electron or 
positron requires that nfld = Q c ±. Thermal dispersion Au/u\ > m±/m p in the upstream 
ions prevents growth for harmonic number n > M s = u±/Au, which cuts the efficiency of 
generation of a non-thermal tail to the pairs' momentum distributions. 
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Simulation parameters 




rrii/m± = 100 Au/u = 0.1 


AVAL 


0.1 


0.2 


Ui/U tot 


0.72 


0.84 


Ui/U_ 


10 


20 


Downstream e ± distribution parameters 




electrons 


positrons 


electrons 


positrons 


f 


0.85 


0.8 


0.6 


0.6 


T' 


20 


20 


19 


19 




120 


110 


120 


100 


N t 


10~ 3 


2 x 10~ 3 


3 x 10- 3 


4 x 10~ 3 


a 


4.5 


4. 


4. 


3.3 


r 


4 


4 


4 


4 


Imax/ Imax 


0.09 


0.1 


0.15 


0.65 


V [%} 


2.5 


4 


4 


4 



Table 3: The table refers to simulations employing a mass-ratio of 100 between the ions and 
the pairs and a finite temperature upstream: Au/u = 0.1 (plots in Fig. 12). In the upper 
part of the table we report the values of the paramaters with which each simulation was 
performed: number ratio between ions and electrons; fraction of the total upstream energy- 
density carried by the ions; ratio between the ions' and electrons' energy density upstream. 
In the lower part of the table, we report the best fit values of the parameters in Eq. 7, 
describing the electrons' and positrons' distribution downstream. In the bottom line we 
report the energy acceleration efficiency for the pairs, estimated (Eq. 8) as the percentage 
of the total upstream energy density U to t converted into non-thermal e ± in each case. 



3. Discussion and Conclusions 



In accord with previous results (Langdon et al. (1988); Gallant et al. (1992)), we find 
full thermalization in pure e~ — e + shocks. The reason is the rapid reabsorption of waves 
at and above the upper hybrid frequency fluH±, which causes the shock to behave as if it 
were collisional - rapid emission of the EM waves by the relativistic cyclotron instability 
at the particle reflection front, with rapid reabsorption of the waves, and exchange of real 
photons, is a mechanism just as capable as the exchange of virtual photons in Coulomb 
collisions in producing local thermodynamic equilibrium. The conclusion, that pure pair 
shocks yield relativistic Maxwellians in the downstream medium, is not limited to the one- 
dimensionality of our model, as is shown elsewhere in 3D simulations by Spitkovsky and 
Arons (in preparation). Those simulations show that at sufficiently low er (er < 1CT 3 is 
low enough), the relativistic Weibel instability (e.g., Yoon and Davidson (1987)) takes over 
from cyclotron instability in the pairs as the deceleration and thermalization mechanism of 
the shock wave, even when the background magnetic field is fully transverse to the flow. 
The physical reasons for this result are discussed in detail elsewhere, in the context of 
the 3D results 1 . However, the essence is that the magnetic fields generated by the Weibel 
instability in an unmagnetized shock saturate in a tangled network of tangled magnetic 
filaments (Milosavljevic & Nakar (2006)), taking a distance kherm ~ 25c/ u p i. These filaments 
scatter and thermalize the incoming particles. Since ltherm/ r Li± ~ 25^/<7i±, Weibel based 
thermalization (a multi-diemensional effect) controls the shock structure for a\± < 0.002. 
For larger magnetizations, the cyclotron effects studied here dominate, somewhat modified 
by mild crinkling of the shock front. 

In the shock wave, the incoming flow continuously injects cold (M s ^> 1) ions into 
the pair plasma, the latter being already heated by the leading shock in the pairs alone. 
These freshly injected cold ions efficiently produce high harmonic waves which are cyclotron 
resonant with the pairs, even though the heated ions that are found further downstream 
from the pair shock do not efficiently radiate high harmonic waves. Magnetic reflection at 
the shock front produces a continuous supply of cold gyrating ions, thus ensuring energy 
transfer from the ions to the non-thermal pairs in the downstream medium, different from 
what is the case in the analogous initial value problem in a uniform medium. 

When the ion density is a small fraction of the pair density, non-thermal acceleration 



1 Other recent simulations of Weibel mediated shocks in e ± and e—p plasmas (Silva et al. (2003); Hededahl 
et al. (2004); Frederiksen et al. (2004); Hededahl & Nishikawa (2005); Jaroschek (2005)) have had insufficient 
size and temporal length (or have been constrained by periodic boundary conditions in the flow direction) 
to show the development of the shock wave; as a result, they have been, in effect, studies of the foreshock. 
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becomes symmetric between electrons and positrons. It also becomes less efficient, leading 
to steeper power-law slopes in the non-thermal part of the spectrum. As shown in Fig. 10, 
the spectra take a Maxwellian form up to energies around E± '-fim±c 2 , then have a power 
law form at higher energies. The energy above which the power law distribution sets in, and 
the slope, depend on Ni/N±. Thus, apparent universality (if it exists) of power law spectra 
in synchrotron sources where relativistic shocks are inferred to play a role in non-thermal 
accleration would be a consequence of universality in the composition of the upstream flow, 
expressed as the ratio Nu/N 1± , if the cyclotron resonance mechanism described here really is 
the dominant non-thermal heating process and if all relativistic shocks are sufficiently close 
to the transverse geometry studied here. 

At present, nothing has been published on the full 3D structure of these shocks 2 . Our 
expectation is that the 3D structure of electron-positron-ion transverse shocks will show a 
transition to unmagnetized physics at a low value of (J^tatai- The value of this transition 
magnetization is quite uncertain, since current neutralization of the ion filaments formed by 
the Weibel instability by electron return currents (Lyubarsky & Eichler (2006)) may greatly 
increase the thermalization length, thus giving more room in a space for cyclotron dynamics 
to ultimately mediate the shock structure. These issues are under investigation. 

Relativistic flows that get their energization from Poynting losses from compact objects 
(possibly all in our list of pulsar winds, microquasar, AGN and GRB jets) may well be in this 
category, since the pair and ion content may both be set by the underlying electrodynamics 
of the compact objects (rotating neutron star, magnetized disk.) Quantification of how close 
shocks have to be to transverse requires consideration of a wider range of shock geometries. 
Investigation of other non-thermal heating mechanisms, such as Diffusive Shock Acceleration 
and magnetic pumping, require higher dimensional simulations. Spitkovsky and Arons (in 
preparation) address both of these issues. 

In the shocks, the relative proportion of electron versus positron acceleration follows 
the ion to pair density, as expected from the changes in wave polarization inferred from the 
linear theory. 

Our study of the shock acceleration process confirms, therefore, the results found by 
Hoshino et al. (1992), extending them to a higher mass ratio. Also confirmed is the ex- 
pectation that signs of electron acceleration would be found if one could consider a plasma 
that contained the same amount of upstream flow energy given to a smaller fraction of ions. 
The overall evolution of the particles is consistent with resonant cyclotron absorption of 



2 Hedcdahl ct al. (2004) in essence study the foreshock of an e — p relativistic shock in 3D with a mass 
ratio m p /m e = 15. 
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magnetosonic waves by the pairs, whose frequencies are harmonics of the ion cyclotron fre- 
quency, being the main process providing non-thermal acceleration. Magnetic reflection of 
the upstream ions in the shock front sets up the unstable ion cyclotron motion which causes 
the wave emission. The non-thermal acceleration efficiencies and power law indices of the 
particle spectra are shown, as a function of the basic governing parameter iVj/JV_, in Table 
2, when the upstream flow is cold. The reduction of the acceleration efficiency for a single 
value of the Mach number M s = 10 <C rrii/m± = 100, with results summarized in Table 3, 
shows the strong effect of large upstream momentum dispersion in the upstream flow. Equiv- 
alently, the absorption of ion cyclotron waves in the shock structure by the pairs is a strong 
non-thermal acceleration mechanism when the upstream flow is quite cold (M s > mj/ra±), 
as results from adiabatic expansion in relativistic outflows from compact objects. 

We conclude by pointing out an interesting feature of this acceleration mechanism. 
The acceleration efficiency and the spectral form of the accelerated particles depend on the 
composition of the flow upstream of the shock. That composition depends on the source 
of the relativistic flow within which the shock wave occurs. Relativistic flows of the type 
considered here are thought to be the consequence of magnetized wind and jet acceleration 
by compact objects, and pair creation within and near such objects. These properties are 
the result of the compact objects' electrodynamics, which constrains some aspects of the 
outflowing plasma composition to particular values - for example, in pulsars and magnetars, 
the ion flux is constrained to be the Goldreich- Julian value. Thus, one might expect a 
correlation between the radiative emission properties of the shocked flow and the emission 
properties of the underlying compact object, as appears to be the case in the X-ray emission 
from pulsar wind nebulae and the spectra of the pulsed X-ray emission from the pulsars 
driving those nebulae (Gotthelf (2003)). 
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Science, and from the forbearance of the taxpayers of California. We accomplished part of 
the work while we enjoyed the hospitality of the Kavli Institute for Theoretical Physics in 
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A. Relativistic Cyclotron Instability of Ring Distributions in a Uniform 

Magnetized e + — e" — p plasma 

We consider an overall neutral e + — e~ — p plasma uniformly distributed in space and 
permeated by a uniform magnetic field. The initial distribution of the plasma in momentum 
space is described as rings, with all species gyrating in the background magnetic field with 
the same Lorentz factor, 7+ = 7_ = j p . For most of the work in this section, we assume 
that all the species are cold, that is, have no dispersion around their gyration momenta. In 
preparation for the shock studies described in §2, we also briefly describe some results of 
the linear instability theory when the ions gyrate with a small momentum dispersion in a 
relativistically hot pair plasma. 

Cold ring distributions are known to be unstable to cyclotron emission and the growth- 
rate of the instability scales proportionally to the particles' cyclotron frequency. Hence, in the 
situation we are presently considering the process of relaxation of the particles' distribution 
functions will show two different stages. 

Initially, the pairs, whose gyration frequencies are much larger than those of the ions, 
will undergo cyclotron instability and tend to thermalize through emission and absorption 
of relativistic cyclotron waves. During this stage, the ions form a static uniform neutralizing 
background: they cannot interact efficiently with the waves produced by the pairs, which are 
too high frequency, and therefore their distribution is expected to be left nearly unchanged 
while the pairs relax. Only at later times the ions will develop the same instability, with 
the delay determined by the difference in the Larmor frequencies of the lighter and heavier 
particles. 

In the case of a plasma with about equal amount of energy carried by light particles 
and magnetic fields (er_ ~ er+ ~ 1), the time-scale for the cyclotron instability in the pairs 
to reach saturation turns out to be of order a few tens of gyration periods, while it is shorter 
for lower magnetization (see Hoshino & Arons (1991)). Hence, if one considers, as we shall 
do in the following, a plasma with magnetization a <, 1, for high enough mass ratios, the 
expectation is that the instability of the pairs will have reached saturation, and they will 
have already almost completely thermalized, before the ion cyclotron instability sets in. 

A.l. Linear theory of the ion cyclotron instability 

The ion distribution function is described by: 

fi(u ± ,u\\) = — !— 5(u ± -u )5(u\\), (Al) 
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where u is the particle four- velocity (u = 7/3, with 7 = 1/ y/l — f3 2 being the particle Lorentz 
factor), uo is the initial value of u and the symbols || and _L are with respect to the magnetic 
field direction. We assume the magnetic field aligned with the z-axis of a cartesian system of 
coordinates B = Be z . Since we are interested in the particles' cyclotron emission the waves 
we are going to consider will have both the electric field and wave vector perpendicular to 
the magnetic field direction (extraordinary mode). If we assume k = ke x the dispersion 
relation can be written as: 

_ _ £ xy£yx ( \o\ 

—•T - £ yy ~ • 

W <^xx 

The general expression for the dielectric tensor in a uniform magnetized plasma can 
be found, for example, following the theory in the book by Krall & Trivelpiece (1973) (see 
also Hoshino & Arons (1991)). We mentioned above that when the ion instability is still 
in its early phases, the pairs will have already exhausted their free energy. This fact leads 
us to expect the fluid approximation to appropriately describe the pairs' contribution to 
the propagation characteristics of the waves whose growth is stimulated by the ions. We 
will therefore describe the pairs' response as that of a relativistically hot fluid with effective 
temperature T± (/cbT ± = P±/N± where P± is the pressure and N± the number density) and 
adiabatic index T. The latter would be T = 4/3 for a 3D relativistic plasma and T = 3/2 
for a 2D relativistic plasma: we confine our attention to the 2D case, since our spatially ID 
simulations do not heat the plasma in the magnetic field direction. 

The relevant components of the dielectric tensor will then be given by: 

n 2 7o f , . d fi (J n (®iU±/lo)) 2 



1 + 27r 2^ 7^T^2 / du \\ du - 



du ± ( v /TT^i/7o - n/ui) 
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where we have used the following definitions: 



f^o B 2 kc uj 

W 



<T > = 7fT L =A 2' a * = ?P' ^ = ?f-' ( A6 ) 

ll* 0s Anm s N s >y Qs c 2 il c0s il c0s 



with 



Q P 0s = \ 4:7rNses = an d £} c0s = 6sB ° = (A7) 

■ irises m s 7 0s c 7 0s 



the plasma and cyclotron frequencies of particles of species s respectively. The functions J n 
that appear in the Eqs. A3-A5 are the ordinary Bessel functions of the first kind of index n 
and argument z = cj3ik/Q ci . 

Moreover, here, 7o« is equal to the initial Lorentz factor of the plasma 70, whereas 7o±, 
entering the quantities that refer to the pairs, is defined by the temperature of the pairs 
after the saturation of their instability: 7 ± = kBT±/m±c 2 . The relation 7 0± < 70 is going 
to hold, since, even if the fraction of energy the pairs lost to waves during the instability 
were negligible, the pressure of a Dirac 5 distribution function centered on 7 corresponds to 
that of a maxwellian with 7 = (7 — l)/2 ~ 7/2. 

If the ions can be treated as a cold fluid when the instability sets in for them, one can 
perform the integrations in Eqs. A3-A5 by parts (see Hoshino & Arons (1991)). 

In general, solutions of the dispersion relation are found numerically and in Fig. 13 we 
plot the results in a few cases relevant for the present study. In all plots we consider a plasma 
with m,i/m± = 2000 and with a baryon fraction Aj/AL = 0.01. The two upper panels refer 
to the case when the ions can be treated as completely cold (analogous to Hoshino & Arons 
(1991)). In the upper left panel the pair plasma is approximated as neutral and the coupling 
between the longitudinal and transverse components of the extraordinary mode is neglected 
e xy = 0. In the upper right panel roots of the complete dispersion equation are shown. 

When the coupling is neglected, at high harmonic numbers, the dispersion relation 
shows two separate unstable branches, corresponding to electromagnetic (cu ~ kc) and mag- 
netosonic waves (u « kc^/ (1/2 + o r ±)/(l + o±))- The growth-rate of magnetosonic waves is 
more or less constant with increasing n while the growth of electromagnetic waves is soon 
depressed at high harmonics. The main effect of the coupling is to further depress the growth 
of electromagnetic waves, whereas the growth-rate of magnetosonic waves stays about the 
same. 



In the lower panels of Fig. 13, we show the solution of the dispersion equation in cases 
when the ion component is not treated as perfectly cold at the onset of the instability. The 
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Fig. 13. — The dispersion relation for waves produced by the ion cyclotron instability. A 
mass ratio of 2000 between the ions and the pairs is used, and the energy balance between 
the different components corresponds to Oi = 0.1 and er_ = 2 ~ er + . Upper left panel: 
the dispersion equation is solved in the approximation that the ions are completely cold and 
that mode coupling is negligible (e xy = 0). Upper right panel: ions treated as cold but 
mode coupling included. Lower left panel: a thermal spread of the ions corresponding 
to Au/uq = 0.01 (see text) is assumed. Lower right panel: a thermal spread of the 
ions corresponding to Au/u = 0.1 (see text) is assumed. While in the upper panels both 
the electromagnetic and magnetosonic wave branches have non-negligible maximum growth- 
rates, only magnetosonic waves are found to grow in the presence of a finite spread in the 
temperature of the ions. 
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plasma composition is the same as above but now a thermal spread of the ion distribution 
is included, with the distribution function in Eq.Al replaced by: 

The effect of the thermal dispersion of the momenta depends on the effective Mach number 
M s = u 0i /Au; we use this terminology, anticipating the occurrence of the ion cyclotron in- 
stability in the magnetosonic shock structures studied in §2, where uoi becomes the upstream 
4-velocity and Aw the upstream dispersion in the ion 4-velocity. We assume M s = 0.01 in 
the left panel, while the case M s = 0.1 is on the right. It is apparent that the growth of the 
waves is suppressed at high harmonic numbers, with the suppression affecting increasingly 
low harmonics the larger the distribution in perpendicular momentum space. 
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